function D = EuDist2(fea_a,fea_b,bSqrt)
% EuDist2 efficiently computes the Euclidean distance matrix
%   D = EuDist(fea_a,fea_b)
%   fea_a:    nSample_a * nFeature
%   fea_b:    nSample_b * nFeature
%   D:      nSample_a * nSample_a
%       or  nSample_a * nSample_b
%   Written by Deng Cai (dengcai@gmail.com)
if ~exist('bSqrt','var')
    bSqrt = 1; 
	% bSqrt=1: Euclidean distance; bSqrt=0: Euclidean distance^2
end
if (~exist('fea_b','var')) || isempty(fea_b)
    aa = sum(fea_a.*fea_a,2);
    ab = fea_a*fea_a';    
    if issparse(aa)
        aa = full(aa);
    end    
    D = bsxfun(@plus,aa,aa') - 2*ab;
    D(D<0) = 0;
    if bSqrt
        D = sqrt(D);
    end
    D = max(D,D');
else
    aa = sum(fea_a.*fea_a,2);
    bb = sum(fea_b.*fea_b,2);
    ab = fea_a*fea_b';
    if issparse(aa)
        aa = full(aa);
        bb = full(bb);
    end
    D = bsxfun(@plus,aa,bb') - 2*ab;
    D(D<0) = 0;
    if bSqrt
        D = sqrt(D);
    end
end